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ABSTRACT 

We have performed two-dimensional hydrodynamic simulations of the collapse 
of isolated axisymmetric clonds condensing via radiative cooling in a primordial 
backgronnd gas. In order to study the development of the so-called “shape- 
instability”, we have considered two types of axisymmetric clonds, oblate and 
prolate clonds of various sizes and with axial ratios of 0.5 < Rc,r/-Rc,z < 2. We 
hnd that the degree of oblateness or prolateness is enhanced dnring the initial 
cooling phase. Bnt it can be reversed later, if the initial contrast in cooling times 
between the clond gas and the background gas is mnch greater than one. In such 
cases an oblate cloud collapses to a structure composed of an onter thin disk 
and a central prolate component. A prolate clond, on the other hand, becomes a 
thin cigar-shape strnctnre with a central dense oblate component. The reversal 
of shape in the central part of the cooled clonds is dne to snpersonic motions 
either along the disk plane in the case of oblate clonds or along the symmetry 
axis in the case of prolate clonds. For a backgronnd gas of = 1.7 x 10®K 
and ri/i = 0.1 cm“^ in a protogalactic halo environment, the mean density of the 
clond gas that has cooled to lO'^K increases to lOOri/i or so, in our simulations 
where noneqnilibrinm cooling is adopted and the backgronnd gas cools too. The 
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spherical Jeans mass of such gas is estimated to be about Mj ~ 5 X 10^ M0. 
In order for cloud mass to exceed the Jeans mass and at the same time in order 
for the thermal instability to operate, the initial cloud size should be around 
1 — 1.5/cooi where /cool is the cooling length. 

Subject headings: galaxy: globular clusters: general - hydrodynamics - instabil¬ 
ities 


1. Introduction 

In many areas of astrophysics, the thermal instability is often invoked to explain the 
condensation of cold dense clouds out of a hot background medium {e.g., Field 1965; Gold¬ 
smith 1970; Defouw 1970; Schwarz et al. 1972; Fall & Rees 1985; Balbus & Soker 1989; 
Vazquez-Semadeni et al. 2000; Koyama & Inutsuka 2002; Kritsuk & Norman 2002). In the 
simplistic picture of the thermal instability, the overdense region surrounded by the hot¬ 
ter background undergoes a “quasi-static compression” in near pressure equilibrium. This 
“near-equilibrium” case of the evolution, however, is valid only when the cloud is small 
enough to adjust to pressure change faster than it cools - its size must be much smaller 
than the distance that sound wave travels in a cooling time, /cool- The collapse of thermally 
unstable clouds in either X-ray cluster cooling flows or protogalactic halos has been studied 
previously by numerical simulations {e.g., David et al. 1988; Brinkman et al. 1990; Kang 
et al. 2000, and references therein). These simulations showed that a spherically symmetric 
cloud cools and undergoes a supersonic compression, if the cloud size is Rc ~ /cool- Such 
supersonic compression leads to a central density increase two to three orders of magnitude 
higher than what is expected from the isobaric compression. On the other hand, a small 
cloud of Rc < /cool cools isobarically and undergoes a quasi-static compression, while a big 
cloud of Rc > /cool cools nearly isochorically with only small density increase. 

While most previous studies on the thermal instability considered the collapse of one¬ 
dimensional (ID) spherically symmetric clouds, Brinkman et al. (1990) simulated the col¬ 
lapse of a non-spherical, elongated blob in the two-dimensional (2D) polar geometry and 
showed the development of a “shape instability”. As the perturbation is compressed by the 
background pressure, the compression wave travels the same distance with the sound speed 
{i.e., I ^ Cg -t) along both the major and minor axes, and the induced infall velocity held is 
not radial. This causes the compressed region to be more elongated and enhances the degree 
of non-sphericity. They showed that an oblate cloud with an initial axial ratio ofb/a = 0.447 
collapses to a hat pancake, and argued that the evolution of oblate clouds becomes similar 
to that of ID plane-parallel collapses. From this simulation one can deduce that the shape 
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instability would also occur in a prolate cloud, resulting in a thin rod-shape condensation. 
In Kang et al. (2000), we studied the effects of different geometries by ID plane-parallel and 
ID spherical symmetric simulations of thermally unstable clouds using a PPM (Piecewise 
Parabolic Method) hydrodynamic code with self-gravity and radiative cooling. We found 
that isotropic compression leads to much higher central density, accompanied by acceler¬ 
ated radiative cooling, in the spherically symmetric case, while the density increases only 
to the isobaric ratio of pdPh = (ThPc/TcPh) in the plane-parallel case. Here p, is the mean 
molecular weight, and the subscripts c and h stand for cool and hot, respectively. Thus, we 
expect that the cloud shape is an important factor in the collapse and evolution of thermally 
unstable clouds. In this contribution, we explore in detail how non-spherical perturbations 
in a primordial environment evolve under the influence of the thermal instability. 

In many models of globular cluster (GC) formation, dense protoglobular cluster clouds 
(PGCCs) are supposed to exist in pressure equilibrium with the hot gas in protogalaxies 
{e.g., Gunn 1980; Brown et al. 1991; Kumai et al. 1993). The model by Fall & Rees (1985), 
which has been most widely adopted to explain the existence of such PGGGs, relies on 
the thermal instability for the formation of PGGGs in a protogalactic halo. In our study 
the parameters have been adjusted so the simulations for the nonlinear development of 
the thermal instability are applicable to the condensation of PGGGs, although the generic 
results should hold to the thermal instability of any objects. We have ignored self-gravity, 
because the gravitational time scale is much longer than the cooling time scale during the 
early collapse stage of PGGGs. But the self-gravity should be important in the gravitational 
fragmentation of PGGGs during the later evolutionary stage when the clouds have cooled 
down to lO^K. 

In next section, we describe our model and numerical method. The simulation results 
are presented in §III, followed by the conclusions in §IV. 


2. HYDRODYNAMIC SIMULATIONS 
2.1. Isolated Clouds in an Uniform Backgronnd Halo 

We expect that inside a protogalactic halo, density perturbations on a wide range of 
length scales exist and flow motions are likely turbulent. While the thermal instability 
under such realistic global pictures can be explored later (c/., Vazquez-Semadeni et al. 
2000; Kritsuk & Norman 2002), here we first take a much simpler and local approach. We 
consider isolated overdense clouds embedded in a hot, uniform, and static background gas of 
Th = 1.7 X 10®K and Uh = 0.1 cm This temperature corresponds to that of an isothermal 
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sphere with circular velocity Vc = 220 km s“^, representing the halo of disk galaxies like 
the Milky Way Galaxy, rih is the background density of hydrogen nuclei. The value of 
n/i = 0.1 cm“^ is chosen as a hducial value, because then spheres of radius Rc ~ /cool would 
have mass scales relevant for GC formation (see §2.3). We assume the ratio of He/H number 
densities is 1/10, so that the gas mass density is given by ph = (2.34 x 10“^^g)n/i. 

The initial density of the overdense clouds, Udoud, is assumed to decrease gradually from 
the center to edge. 


^cloud 


nh 


1 + 5 cos 




( 1 ) 


for + (z/Rc^z)^ < 1, where i?c,R and i?c,z are the cloud radii along the R and 

2 : axes, respectively. The ratio of i?c,R/-Rc,z determines the shape of the initial clouds 
(see §2.5). The initial temperature throughout the clouds is set by the isobaric condition, 
i.e., Tcioud = Th{nh/ncioud)- The amplitude of initial density perturbations, 5, is a free pa¬ 
rameter that determines the density contrast between the cloud center and the background. 
If this amplitude is linear {i.e., 5 +; 1) and if there are no heat sources that balance the 
radiative cooling of the background gas, then both the cloud and the background gas cool 
together and the thermal instability would not have enough time to grow. Because the 
cooling time scales as fcooi oc (1 + 5)“^, the contrast in fcooi between the cloud and the back¬ 
ground is small for small 6. In real protogalaxies, however, the halo gas would be heated 
by possible energy sources such as supernova explosions, shocks, and etc. If the background 
gas maintains a high temperature owing to those energy inputs, perturbations would grow 
approximately under the isobaric condition until they become nonlinear. Also these heating 
processes likely induce turbulent flow motions and non-linear density fluctuations in the halo 
medium as well. In numerical simulations the overall evolution proceeds faster for larger val¬ 
ues of S during the linear phase, but it becomes almost independent of the initial values of 
5, once perturbations become non-linear. Since we do not include any background heating 
processes in our simulations, we need to have a reasonably large density contrast in order to 
see the nonlinear growth and to expedite the simulations. Thus, we begin with 5 = 1 for all 
models, which in fact would be consistent with turbulent nature of the halo medium. 


2.2. Radiative Cooling Rates 

The key idea of the GG formation model based on the thermal instability, originally 
suggested by Fall & Rees (1985), is that the characteristic mass scale of GGs, Me ~ 10® Mq, 
can be explained by the imprinting of the Jeans mass of the gas clouds at T = lO^K that 
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have cooled from a hot halo gas in pressure equilibrium. If the clouds were allowed to cool 
well below 10*^ K in a time scale shorter than the free-fall time scale, they would not retain 
the memory of the imprinted mass scale. So in this work we consider the cases where the 
radiative cooling is ineffective below lO^K. If the halo gas had been enriched by the metals 
from the first objects, or if H 2 molecules have formed efficiently via gas phase reactions, then 
the gas would have cooled well below 10^ K before the Jeans mass was imprinted (Shapiro 
& Kang 1987). Thus, we consider a primordial gas with H and He only, and we assume 
that the formation of H 2 molecules is prohibited due to UV radiation from central AGNs or 
diffuse background radiation (Kang et al. 1990), resulting in zero cooling for T < 10"^K. 

Here, we dehne the cooling rate as 

^{T,nH)=L{T)nl, ( 2 ) 

where A is the energy loss rate per unit volume and uh is the number density of hydrogen 
nuclei. In general, the cooling rate coefficient, L{T), is a function of temperature as well 
as ionization fractions that can be dependent on the thermal and ionization history of gas. 
When a hot gas cools from T > 10®A', in particular, it recombines out of ionization equilib¬ 
rium, because the cooling time scale is shorter than the recombination time scale (Shapiro & 
Kang 1987). So in order to calculate the non-equilibrium cooling rate accurately, the time- 
dependent equations for ionization fractions should be solved, which can be computationally 
expensive. Fortunately, however, the non-equilibrium cooling rate for the gas cooling under 
the isobaric condition becomes a function of temperature only, if the initial temperature is 
high enough to ensure the initial ionization equilibrium {e.g., T > 10®K) and if only two- 
body collisional processes are included. In that case, L{T) can be represented by a tabulated 
form as a function of temperature only. We adopt, as our standard cooling model, the non¬ 
equilibrium radiative cooling rate for a zero-metalicity, optically thin gas that is calculated 
by following the non-equilibrium collisional ionization of the gas cooling from 10^ ®K to 10^ 
K under the isobaric condition (Sutherland & Dopita 1993). We set L{T) = 0 for T < lO^K 
and, in addition, set the minimum temperature at Tmm = lO^K. 

In order to explore the effects of different cooling rates, we have also calculated several 
models with the following cooling rates, in addition to the standard cooling model (NEQmO): 
1) the CIEmO model: the collisional ionization equilibrium cooling for a zero-metalicity 
gas, 2) the NEQml model: the non-equilibrium cooling rate for a gas with the metalicity, 
Z = 0.1 Zq. The NEQml model is not consistent with our assumption of zero cooling rate 
below T = 10^, but has been included for comparison. Figure 1 shows the cooling rate 
coefficients for these cooling models. 

As the central part of the cloud cools, a steep temperature gradient develops between 
the cloud and hot background medium and the thermal conduction can become operative 
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there. Brinkman et al. (1990) showed, however, that their simulation results differ by <3% 
when the reduced thermal conductivity (Gray & Kilkenny 1980) was included. Thus we 
ignore the thermal conduction in our calculations, although possibly it may be important. 


2.3. Cooling Time and Length 


We dehne the cooling time as 


Gool 


3/2ntotkT 

A 


(3) 


where e = p/{'y — 1) is the internal energy per unit volume, and Utot is the number density 
of ions and electrons. With the standard NEQmO cooling model, the cooling time for the 
background halo gas is 

Gi,h = 2.0 X 10^ yrs ■ { 3 )"^ (4) 

0.1 cm 

with Th = 1.7 X 10®K. Figure 1 shows the cooling time for a gas cooling under the isobaric 
condition {i.e., ntotT= constant) for three cooling models. The cooling time for the gas at 
the cloud center is given by 

^ _ tcl,h Gl,h 

- (T^ = — 

with ndoud = nh{l + S), Tdoud = Th/{1 + 5) and 5 = 1. We note the cloud center cools in a 
time scale of ~ 0.5tci,h, while the background halo cools in ~ 2tci,h- 


We dehne the cooling length as the distance over which the sound wave of the hot halo 
gas travels within one cooling time of the perturbed gas, 

/cool = Ch ■ tci,c = l.OSkpc ■ —r)"\ (6) 

U.l cm 

where Ch = 198 km s“^ is the sound speed of the hot gas of T/^ = 1.7 x 10®K. The dynamics 
of radiatively cooling clouds are characterized by the cloud size relative to the cooling length 
(Fall & Rees 1985; David et al. 1988; Kang et al. 2000). 


The mass contained within a cloud of radius, Rc = /cool, is 

Af„„i = = 1.75 X 10 ' Ms ■ ( "'■ 3 )-'. (7) 

3 0.1 cm 

For Uh = 0.1 cm“^, the clouds with the size Rc ~ (0.4 — l)/cooi have the mass scales of 
Tfcioud ~ 10® — 10 ’^ M 0 that are relevant for PGCCs. Because of the n'^‘^ dependence, for 
the background halo density that is much larger or much smaller than our hducial value, 
the characteristic cooling mass would be too small or too big, respectively, for the cooled 
perturbations to become PGGGs. 
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2.4. Numerical Method 


The gasdynamical equations for an axisymmetric system in the cylindrical coordinate 
system, {R,z,(f)), including radiative cooling are written as 


where e = 
u = u^ + u 


dp Id d 

dipun) , 1 (9 ^ ^ 2 ^ , <9 , ^ dP ^ 

-IT + M 


d{pu^) 1 d 
dt ^ RdP 
d{pe) 1 d 


d 


dP 


(Rpunu,) + 7r(P«J - 


dz 

d 


dz 


dt 


+ P)] + -^Wzipe + P)] = -A, 


ejf> + (l/2)'((^ is the total energy of the gas per unit mass, R 
and the rest of the variables have their usual meanings. 


( 8 ) 

(9) 

( 10 ) 

( 11 ) 

\/x^ + y^, 


To solve the hydrodynamic part, we have used an Eulerian, grid-based hydrodynamics 
code based on the “Total Variation Diminishing (TVD)” scheme (Ryu et al. 1993). The cylin¬ 
drical geometry version has been used for 2D simulations. For ID comparison simulations, 
the spherical geometry version has been used. The TVD scheme solves a hyperbolic system 
of gasdynamical conservation equations with a second-order accuracy. Multidimensionality 
is handled by the Strang-type dimensional splitting (Strang 1968). 


After completion of the hydrodynamic part updating hydrodynamical quantities from 
the time step P to radiative cooling is applied to the thermal energy as a separate part. 
If we were to update the thermal energy density by the following explicit scheme, 

= 6” - AAf = 6”(1 - ^), (12) 

^cool 

the time step size should be smaller than the cooling time, that is. At < ■ Here 

is estimated from the time averaged hydrodynamic quantities, (g" -|- g”’''^)/2. This 
explicit integration scheme would be extremely expensive, when the cooling time scale is 
much shorter than the hydrodynamical time scale. For that reason, we rewrite the cooling 
part of the thermal energy equation as 


(iln(e) 


A 

e 


dt 6 Pool 

In our numerical code we integrate this equation as 

Ath , 


(13) 


gn+l _ gn . 


,n+l/2f' 

^cool 


(14) 
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assuming that tcooi is constant over the hydrodynamic time scale Ath {e.g., LeVeque 1997). 
Strictly speaking, Equation (14) converges to Equation (12) only when Ath is much smaller 
than tcooi- When Ath ^ however, the gas quickly looses most of the thermal energy 

via radiation and the temperature approaches to the specihed minimum value (Tmin) during 
one hydrodynamic time step, if we were to integrate Equation (12) with a small timestep 
size of tcooi for many steps. But this behavior is emulated reasonably well by integrating 
Equation (14) with Ath for one step, which effectively lowers the gas temperature to T^in- So 
the end results of both integration schemes would be qualitively similar, although could be 
quantitatively somewhat different. Especially at different spatial grid resolutions, the detail 
thermal history of the gas could be different during the cooling phase, but once cooled the 
hnal structure should be roughly similar. With this integration scheme, gas cooling can be 
followed, although approximately, with the hydrodynamic time steps, even when the cooling 
time scale is shorter than the hydrodynamic time scale. 

Our simulations start at t = 0 with the clouds at rest in pressure equilibrium (Phaio = 
Pcioud), and cease at t = 2fcooi when the background gas has cooled to lO^K. The standard 
mirror condition has been used for the reflecting boundaries at i? = 0 and z = 0, while it 
has been assumed that flows are continuous across the outer boundaries. We have used only 
one quadrant of the cylindrical coordinate system. 

In our simulations the clouds cool and collapse due to the compression by the background 
pressure to the size that is 1/10 or so of the initial cloud radii. In order to study the detail 
structure of the collapsed clouds, we have devised a special grid that consists of an inner fine 
zone with uniform grid spacing and an outer coarse zone with expanding grid spacing. The 
inner zone has 1000^ uniform cells with AR = Az = Lfine/1000, where Lfine = 1-6/cooi- The 
outer zone has 1240^ — 1000^ cells covering the rest of the simulated region outside the inner 
hne zone, with the cell spacing that increases outwards as AR^/AR^~^ = = 1.05 

for i = 1001, ...1240. With the expanding grid the outer boundaries are located far away 
from the central cloud where most activity occurs. 

The physical variables are expressed in units of the following normalization both in 
the numerical code and in the plots presented below: to = Pi,h = 2.0 x 10^ years; rp = 
/cool = 1.05kpc; Mo = /cooi//ci,h = 50.8kms“^; po = 2.34 x lO^^^g ■ Uh] Po = Po^l = 6.03 x 
10“^^erg cm“^. 
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2.5. Shape Parameter 


In order to quantify how the cloud shape evolves, we dehne the “shape parameter” as 


S 


R 


e.R 


R. 


e,z 


(15) 


where i?e,R and i?e,z are the “effective” radii along the R and ; 2 -axes, respectively. Here, 
the “effective” radius is dehned as the radius where the gas density decreases to a half of 
the central density, that is, p{Re) = (l/2)pc- The initial values of the shape parameter for 
models considered are = 1 for spherical clouds. Si = 5/6 and 1/2 for prolate clouds, and 
Si = 6/5 and 2 for oblate clouds. The second column of Table 1 shows the initial values of 
the shape parameter for each model. 


3. SIMULATION RESULTS 
3.1. Spherical Symmetric Calculations 

We have first calculated spherical symmetric collapses both with ID and 2D codes, in 
order to show their behavior and to compare it with that of non-spherical collapses, as well 
as to test the performance of the 2D code. According to our previous ID simulations which 
employed a different code (the PPM code) (Kang et al. 2000), the evolution of spherical 
clouds collapsing via the thermal instability can be classified by the cloud size as follows: 1) 
{Rc/lcooi) < 0.2, the isobaric compression regime, 2) 0.2 < {Rc/lcooi) ^ 1 — 1.5, the supersonic 
compression regime, and 3) {Rc/lcoo\) ^ 1 — 1-5, the isochoric cooling regime. 

The following four different sizes have been considered in our spherical symmetric calcu¬ 
lations: three supersonic cloud models, Sll with Rc = 0.4/cooi, Mil with Rc = 0.8/cooi, and 
Lll with Rc = 1.6/cooi, one isochoric cloud model, XI 1 with Rc = 3.2/cooi (see Table 1). The 
largest cloud model, Xll, has a cloud mass too large for PGCCs, but it has been included 
for comparison. We hrst note that the results from the current ID simulations are consistent 
with those from our previous simulations in Kang et al. (2000), although different numerical 
codes are used. While the CIEmO cooling model was adopted in Kang et al. (2000), here we 
have adopted the standard NEQmO cooling model which has lower peaks around both the 
H and He Ly a line emissions (see Figure 1). As a result, gas cools less rapidly and central 
density increases to lower values in the current simulations. 

The results from 2D simulations of spherical collapses are mostly similar to those from 
ID simulations. Figure 2 shows the distributions of the gas density, the velocities in the R 
and directions and the pressure along the diagonal line of R = z from the 2D simulation 
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for the Lll model with Rc = 1.6/cooi- The distributions are sampled at t/tci,h = 0.2 (solid 
line), 0.6 (dotted line), 1.0 (dashed line), 1.4 (long dashed line) and 1.8 (dot-dashed line). As 
the central gas cools and loses pressure, the high pressure background compresses the cloud 
and induces an infall flow for t/tci,h ^ 0.6. The flow velocity increases up to {u/ch) ~ 0.36, 
which is much larger than the sound speed of the cooled gas, [cdCh) ~ 0.08. So the accretion 
flow is supersonic. After t = 0.6tci,h, the infall flow bounces back as an accretion shock and 
the shock moves out. Then, the cloud expands outward slowly due to the high pressure 
at the central region. One can see that the central density increases up to pdPo ~ 10® 
just before the shock bounces back, and then decreases to pd Po ~ 10^'® as the collapsed 
cloud expands. The shape of the collapsed clouds from 2D simulations becomes slightly 
rectangular, especially for the models with small sizes (see Figure 8 below). The deviation 
from the spherical symmetry, measured with the shape parameter (S' — 1), ranges from less 
than 1% (in the Xll model) to ~ 20% (in the Sll model). 


3.2. 2D Axisymmetric Calculations 

For the non-spherical, axisymmetric clouds, four different sizes are considered, as for the 
spherical clouds: i?c,max = 0.4/cooi for small clouds, i?c,max = 0.8/cooi for medium size clouds, 
-Rc,max = T6/cooi for large clouds, and i?c,max = 3.2/cooi for very large clouds. Here, i?c,max = 
max[i?c,R, i?c,z]- The initial parameters of the model clouds considered are summarized in 
Table 1. 

Figures 3 and 4 show the time evolution of the density distribution of the large prolate 
cloud, the L56 model with S\ = 5/6, and the large oblate cloud, the L65 model with S\ = 6/5, 
respectively. Figure 5 shows the flow velocity held superimposed on the density contour maps 
of the L56 (left panels) and L65 (right panels) models at t = 0.8tci,h (top panels) and at 
t = 1.6tci,h (bottom panels). As gas cools catastrophically, the clouds implode near the 
center and supersonic anisotropic accretion how is induced. The compression continues until 
t ~ 0.8tci,h when infall how is rehected at the center and an accretion shock forms. This 
“shock formation time” depends on the initial cloud size, as the infall how turns around later 
in larger clouds. At the shock formation time, the infall how held is strongly anisotropic, 
preferentially parallel to the z = 0 plane for the prolate cloud and parallel to the 2 :-axis for 
the oblate cloud. As a result, the degree of prolateness or oblateness is enhanced. So the 
prolate cloud collapses to a very thin rod shape, while the oblate cloud collapses to a hat 
disk, at t ~ 0.8tci,h in the L56 and L65 models. The central density peaks at the time of 
shock formation, and then slowly decreases afterwards. After that time, the accretion shock 
halts the infall how, as shown in the bottom panels of Figure 5. Inside the accretion shock. 
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there exist residual infall motions with the preferred direction different from that of early 
accretion before t = 0.8tci,h, that is, along the 2 ;-axis for the prolate model and along the 
z = 0 plane for the oblate model. So now the contraction proceeds along the ; 2 -axis toward 
the center, which in turn generates the radial outflow along the z = 0 plane inside the cooled 
cloud in the prolate model. In the oblate model, on the other hand, the contraction occurs 
along the z = 0 plane toward the center and induces the bipolar outflow along the z-axis. 
This “secondary” contraction leads to a central bulge that has a shape reverse to the initial 
shape of the cloud. Thus, the prolate cloud results in a structure that consists of a cigar¬ 
shaped outer component and a central oblate component. On the other hand, the oblate 
cloud collapsed to a structure of a pancake-shaped outer component and a central prolate 
component. 

In the 2D simulation of Brinkman et al. (1990), they considered an oblate cloud (S', ~ 
2.2) with small density contrast, 5 = 0.1, in a background with Th = lO^K and nt = 0.1cm“^, 
hotter than ours. They showed that the cloud collapses to a flat structure with density 
enhancement of ~ 100 at t ~ 3 in units of their cooling time. Because their simulation 
started with a small 5, it ended just after the formation of a nonlinear pancake-like structure 
that corresponds to the structure at t ~ 0.9tci,h in our simulation (see Figure 4). Asa result, 
they did not see the development of the inner central prolate condensation. Whether collapses 
end at the cigar/pancake formation stage or they continue to form the inner oblate/prolate 
components depends mainly on the cloud size ratio, Rc/lcoo\, and the contrast of cooling 
times, tci,c/Ai,h- Thus, these two parameters along the initial shape parameter. S',, determine 
the dynamics and mass distribution of the collapsed clouds. 

Figures 6 and 7 show the evolutionary sequences for the medium size prolate cloud, 
the M12 model with Si = 1/2, and the medium size oblate cloud, the M21 model with 
Si = 2, respectively. Due to a higher degree of initial non-sphericity, these models display 
much stronger shape instability than the L56 and L65 models. The shock formation time is 
t ~ 0.6tci,h, slightly earlier than that for the L56 and L65 models because of smaller sizes. 
The inner components due to the secondary contraction grow less signihcantly, compared to 
those of the L56 and L65 models. 

To illustrate how the initial shape parameter and cloud size affect the evolution, we 
plot in Figure 8 the density contour maps at t = 1.6tci,h for the models among listed 
in Table 1 which adopt the standard cooling. For all models, the regions bounded by 
[—0.5i?c,R, 0.5i?c,R] X [—0.5i?c,z, 0.5i?c,z] are shown. So for example, [-0.4, 0.4] x [-0.4, 0.4] in 
units of /cool is shown for the M models. Except in the very large clouds (the X models), the 
shape reversal is observed in the central part of cooled clouds; that is, the formation of the 
inner oblate (prolate) component in the prolate (oblate) models. Spherical models are shown 
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to demonstrate how well the spherical symmetry is conserved in 2D simulations (§3.1). It 
is interesting to see that the collapses of the clouds with initially even 20% deviation from 
the spherical symmetry (i.e., M56, M65, L56 and L65) are very different from the spherical 
collapses. For a given value of Si, the models with smaller cloud sizes have more signihcant 
inner components, since the shock formation occurs earlier. For a given cloud size, on the 
other hand, the models with S'ds closer to one form more significant inner components. In 
the very large cloud models, the sound crossing time along the long axis (i?c,max = 3.2/cooi) 
is much longer than the cooling time. Hence, the spherical cloud (the XI1 model) cools 
and collapses without forming an accretion shock. However, the X12 and X21 models have 
the minor-axis radius, i?c,min = 1-6/cooh so they are able to collapse supersonically along the 
minor-axis forming an accretion shock at f ~ fci,h- As a result, hy t = 1.6tci,h the prolate 
cloud (the X12 model) has been developed into a rod with the central density pdPh ~ 10^-'^, 
while the oblate cloud (the X21 model) into a pancake with pdPh ~ 10^'^. As in the sim¬ 
ulation by Brinkman et al. (1990), the inner central components do not display the shape 
reversal by the time t ~ 1.6fci,h in these models. If we adopt a higher initial density contrast 
or keep the background gas at constant pressure, however, even the X models could have 
developed the inner central components with the reversed shape. 

In order to look at the density enhancement in a quantitative way, we present in Figure 
9 the density line prohles of the collapsed clouds at t/tci,h =1.8 for the initially prolate (left 
panels) and oblate (right panels) clouds. The dotted (dashed) lines show the prohles along 
the i?-axis (; 2 -axis), while the solid lines show the prohles of the 2D spherical models with 
the same size {i.e., the Sll-Xll models) along the diagonal line oi R = z a. function 
of the radial distance r = dR^ + z^. In the prolate models the dashed lines reveal the 
thin rod-shape component along the 2 :-axis, while the dotted lines show the inner oblate 
component along the i?-axis, except in the X12 model where the inner oblate component 
does form. In the oblate models, on the other hand, the dotted lines show the flat pancake- 
shape component along the i?-axis, while the dashed lines reveal the inner prolate component 
along the ; 2 -axis, except in the X21 model. Once again the shape reversal is observed in the 
central region (r < 0.1 —0.2/cooi) of the cooled clouds, except in the X models. The outer rod- 
shape component in the prolate models and the outer pancake-shape component in the oblate 
models have the density close to the isobaric ratio, pcooi/ph ~ {ThPc/TcPh) ~ 10^'®. The inner 
oblate/prolate components have the mean density somewhat higher {pcooi/Ph ~ in 

the M and L models. As mentioned before, the central density peaks at the time of shock 
formation and then decreases afterwards as the clouds expand, so the density distribution 
changes in time. If we compare the central density of different models at a given time after 
the inner shape reversal appears, the oblate models have higher values than the prolate 
models, except in the X models. 
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Figure 10 shows the evolution of the shape parameter, S, for the models with Si = 1/2 
(prolate) and 2 (oblate). During the hrst “shape instability” stage, the shape parameter 
is determined by the cigar-shape component or the pancake-shape component. During the 
“secondary infall” stage after the shock formation, however, it represents the shape of the 
inner components in most models (except in the S12 model). In the oblate clouds (right 
panels), as the clouds collapse to flat pancakes, the value of S increases to ^ 1 nntil the 
shock formation time. Snch time when S reaches the maximum values scales as the sound 
crossing time (4c = -Rc,min/ch), which is proportional to the cloud size. Afterwards, the 
value of S decreases, and it becomes smaller than one due to the formation of the inner 
prolate components, except in the model X21. In the prolate clonds (left panels), the value 
of S decreases to 1 as the clouds collapse to thin rod shapes. After the shock formation 
the value of S increases, and it becomes greater than one in the M12 and L12 models as 
the inner oblate components grow. In the small S12 model, the inner oblate component 
and the onter cigar-shaped component have the similar density, so the effective radius along 
the 2 :-axis represents the length of the cigar component rather than the radius of the inner 
oblate component, resulting in S' = Re,R/Re,z < 1- For this model, we have also calcnlated the 
second shape parameter S'loo, which is dehned as the ratio of the radii where p(i?ioo) = lOOp/^. 
As shown in Figure 10, this second shape parameter becomes S'loo > 1 at t > 1.6tci,h, 
indicating that the inner oblate component indeed forms in the S12 model. In the X models, 
S remains either less than one or greater than one, as expected. 


3.3. Cloud Mass and Mean Density 

In the Fall and Rees model for the formation of PGCCs where the thermal instability 
is assumed to proceed quasi-statically, the “critical mass” dehned for an isothermal sphere 
conhned by an external pressnre (McCrea 1957) 

M,, = 1.18(^] G-^l‘^pl^l‘^ (16) 

\iamH) 

was adopted as the minimum mass for gravitationally unstable clonds. For Tc = 10^ K and 
Ph = 3.28 X 10“^^dyne cm“^, M^r = 2.8 x 10® Mq. However, this simple picture should be 
modihed for the following reasons: 1) the collapse is not qnasi-static and the infall hows can 
become snpersonic, so the compressed clonds are bonnd by the ram pressnre of the infall 
hows rather than the backgronnd pressure, 2) the cooled compressed clouds (PGCCs) may 
have turbulent velocity helds (see Figure 5), 3) the mass distribution of PGCCs can be very 
complex, rather than spherical (see Figures 3-8). 

As an ehort to obtain a better estimation on the PGCC mass, we have estimated the 
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Jeans mass of the cooled gas in our simulations as follows. First, we have calculated the 
mass of the gas with T = lO^K, M 4 . The right panels of Figure 11 show how M 4 increases 
with time in our models. Most of the cloud gas cools to 10"^K by t = 1.4tci,h, so by this 
time M 4 becomes comparable to the initial cloud mass, Mdoud- We present the values of M 4 
at t = 1.4tci,h in the fifth column of Table 1. Next, from the mass and volume of the cells 
with T = 10"^K, we have calculated the mean density of the central region, < >. The 

values of < p 4 > are shown in the left panels of Figure 11. Although the central density 
can increase to values much higher than the isobaric compression ratio of lO^'®, the mean 
density, p 4 , increases only up to that ratio and then decreases, as the background gas cools 
and the cooled cloud expands. The values of < p 4 > at t = 1.4tci,h are listed in the last 
column of Table 1. 

Finally, although it may not be a good approximation either because of the reasons listed 
above, we have computed the Jeans mass of an isothermal uniform sphere with temperature 
Tc and density pc as (Spitzer 1979) 

A/, = pA5 = 5.46(^) ' (17) 

by adopting = 10"^K and Pc =< Pa >■ The values of such Jean masses are plotted in the 
right panels of Figure 11, and listed in the 8 th column in Table 1. In the spherical model, 
the cooled gas in the Lll cloud has mass greater than the Jeans mass, i.e., M 4 > Mj. In 
the non-spherical models, the total mass of the cooled gas in the L56, L65, and L21 clouds 
exceeds the Jeans mass. These clouds might become gravitationally unstable after cooling 
to lO^K. In those models, the clouds might imprint the Jeans Mass of Mj ~ 5 x lO'^ Mq. 
Assuming a star formation efficiency of order 10 %, this mass estimate is a bit too large for 
the characteristic mass scale of the GC mass distribution. This mass scale, however, is much 
larger than the previously estimation by Kang et al. (2000). It is because the cloud density 
increases only to pc ~ lOOp/i in the current simulations where the slower nonequilibrium 
cooling rate has been adopted. In any case, our estimation for the gravitationally unstable 
mass scale should be taken to be very rough, and would depend on the model parameters 
including the density and temperature of the background. In addition, self-gravity has been 
ignored in our simulations, which is expected to be important in the later stage of the cloud 
evolution, as pointed in §1. Because the Xll, X12 and X21 models either have small density 
enhancement or have extremely non-spherical mass distribution, we argue that it would not 
be useful to make the comparison between M 4 and Mj for those models. 
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3.4. Different Cooling Model 

Finally, we have calculated the L56 and L65 models with different coolings, to see their 
effects. They are labeled as EL56 and EL65 for the models with the CIEmO cooling and 
as ML56 and ML65 for the models with the NEQml cooling (see also Table 1). Figure 12 
shows the density contour maps at the shock formation time {t = 0.8tci,h) and at a later 
time {t = 1.4tci,h)- In the ML56 and ML65 models where cooling is enhanced by about a 
factor of 3.5 at T ~ 10®K due to metals, the cloud gas has cooled to T = 10"^ K already at 
t = 0.2tci,h and the compression wave steepens into a pair of shocks (reverse and forward) 
hj t = 0.8tci,h- This produces an elongated ring-like structure in the 2D image, which is in 
fact a dense elongated shell in 3D. This structure along with the pair of shocks collapses at 
the center and then an accretion shock bounces back before t = tci,h- The gas of the EL56 
and EL65 models with the CIEmO cooling cools faster, and so we see the resulting structure 
is more compact than that in the models with the standard NEQmO cooling. As a result, 
the density enhancement reaches < >= 360 — S60ph hy t = 1.4tci,h, which are greater 

than those found in the models with the standard cooling. So the spherical Jeans mass of 
the cooled clouds becomes ~ 10'^ Mq, which is more consistent with the characteristic mass 
scale of GCs (see Table 1). 


4. SUMMARY 

In many astrophysical problems, a two-phase medium may develop by the formation 
of cold dense clouds via the thermal instability. In order to explore how the non-spherical 
shape affects the collapse of thermally unstable clouds, we have performed 2D hydrodynam- 
ical simulations in the cylindrical geometry with the radiative cooling rate for a primordial 
gas. Although we have selected for our simulations a physical environment relevant to a 
protogalactic halo, the overall simulation results can be applied to the collapse of thermally 
unstable clouds of any physical scales, as long as the cooling curve has the pertinent char¬ 
acteristics for the thermal instability. 

The collapse of non-spherical clouds depends mainly on the following factors. 

1) The ratio of the cloud size to the cooling length, Rc/lcooi- As shown in the previous ID 
spherical simulations {e.g., David et al. 1988; Brinkman et al. 1990; Kang et al. 2000), clouds 
of Rc ~ /cool undergo a supersonic compression, resulting in high density enhancements. 
Small clouds with Rc < /cool cool nearly isobarically through a quasi-static compression, 
while large clouds with Rc > /cool cool nearly isochorically. Hence, we have focused on the 
clouds with Rc ~ /cool- 

2) The initial density contrast between the cloud and the background, (1 -f J): Note that 
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the contrast in cooling times, tc\,c/tc\,h, is initially proportional to (1 + 6)~‘^ for an isobaric 
perturbation. Since we have not included any heating source that maintains the temperature 
of the background gas, the growth of nonlinear structures depends on the value of 6. For 
example, with 5 -C 1, only the first stage of nonlinear growth can develop before the back¬ 
ground gas itself cools down. With 5 ~ 1, however, more complex structures form after the 
initial emergence of the shape instabilities. So we have started our simulations with <5 = 1. 
3) The degree of non-sphericity, S = Rc,r/Rc,z- It is obvious that the degree of deviation 
from the spherical symmetry determines the dynamics of the infall flows and the mass distri¬ 
bution of the collapsed clouds. So we have considered both the prolate clouds with Si = 1/2 
and 5/6 and the oblate clouds with S', = 6/5 and 2, in addition to the spherical clouds 
(^i = !)• 

The collapse in our simulations can be described by two distinct stages: the hrst shape 
instability stage during which the non-sphericity grows due to the initial infall, and the sec¬ 
ondary contraetion stage during which the infall occurs predominantly along the direction 
perpendicular to the initial infall flows and a “shape reversal” occurs. Even with initially 
only 20 % deviation from the spherical shape {i.e., Si = 5/6 or 6/5), a strong shape insta¬ 
bility occurs, so the prolate clouds are compressed to thin rods and the oblate clouds are 
compressed to flat pancakes due to strongly anisotropic infall flows. The degree of prolate¬ 
ness or oblateness, however, is enhanced only during the initial shape instability phase up to 
the formation of accretion shocks. Afterwards, secondary infall motions are induced, domi¬ 
nantly along the 2 :-axis for the prolate clouds and along the z = 0 plane for the oblate clouds. 
This secondary contraction parallel or perpendicular to the ; 2 -axis induces, within the cooled 
clouds, the radial outflows in the prolate models or the bipolar outflows in the oblate models, 
resulting in the inner central bulges with the mass distribution opposite to the initial shape. 
As a result, initially prolate clouds collapse to a system that consists of an outer cigar-shaped 
component and a central oblate component, while initially oblate clouds collapse to a system 
that consists of an outer pancake-shaped component and a central prolate component. 

The central density of the collapsed clouds increases until accretion shocks form at the 
end of the first shape instability stage. And then it gradually decreases as the clouds expand, 
since the central pressure is higher than the background pressure. The central density in our 
simulations has turned out to be much lower than that in the previous simulations. In the 
secondary contraction stage, the mean density of the outer pancake or thin-rod components is 
similar to the background density times the isobaric compression ratio of {ThHc/Tcp,h) = 10^'^, 
while that of the inner components reaches only up to an order of magnitude higher than 
that. We note that the density enhancement depends on the radiative cooling rate. It would 
be higher in the simulations with larger cooling rates. We have adopted as the standard 
cooling model, the cooling rate of a primordial gas based on the non-equilibrium ionization 
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fraction tabulated by Sutherland & Dopita (1993). This cooling model has lower rates near 
H and He Ly a line emission peaks than the cooling model based on the collisional ionization 
equilibrium which was adopted in our previous ID simulations (Kang et al. 2000). So we 
have found smaller density enhancements in the current simulations. 

For the protogalactic halo environment considered here, ph = 2.34 x 10“^®g ■ cm“^ and 
Th = 1.7 X lO^K, the spherical Jeans mass of the cooled clouds is about Mj ~ 5 x lO'^ Mq. 
This mass scale is somewhat large for the mass of PGCCs which fragment to form GCs. 
But this is based on a very rough estimation which would depend on model parameters. In 
addition, in a realistic halo environment, the halo gas may be heated by the stellar winds, 
supernova explosions, shock waves and etc. If the halo can maintain the high temperature 
and continue to compressed the PGGGs, the density of PGGGs would have increased more, 
resulting in a smaller Jeans mass. Finally, we have considered here a static halo of uniform 
density. However, protogalactic halos are likely clumpy and turbulent, and the PGGGs may 
have formed in such environment. We leave all these issues, along with extension into 3D, 
to future works. 
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Table 1. Initial Parameters for Model Clouds 


Model 


Cool 

-Rc,z/^cool 

-Rc,max 

d^cloud 

M 4 

Mj ^ 

< P4>‘^ 





(kpc) 

(106 Mo) 

(106 Mo) 

(106 Mo) 

(2.34 X 10-24g . Uh) 

Sll 

1.0 

0.4 

0.4 

0.421 

1.489 

1.156 

45.222 

96.125 

S12 

0.5 

0.2 

0.4 

0.421 

0.376 

0.455 

47.744 

86.236 

S21 

2.0 

0.4 

0.2 

0.421 

0.745 

0.856 

49.797 

79.271 

Mil 

1.0 

0.8 

0.8 

0.842 

11.858 

8.940 

32.831 

182.367 

M12 

0.5 

0.4 

0.8 

0.842 

2.979 

2.687 

47.934 

85.557 

M56 

5/6 

0.67 

0.8 

0.842 

8.243 

8.059 

65.272 

46.139 

M65 

6/5 

0.8 

0.67 

0.842 

9.881 

7.935 

65.407 

45.949 

M21 

2.0 

0.8 

0.4 

0.842 

5.929 

5.928 

50.757 

76.304 

Ill 

1.0 

1.6 

1.6 

1.68 

94.622 

71.112 

58.688 

57.074 

L12 

0.5 

0.8 

1.6 

1.68 

23.715 

19.752 

40.654 

118.941 

L56 

5/6 

1.33 

1.6 

1.68 

56.741 

55.689 

39.360 

126.890 

L65 

6/5 

1.6 

1.33 

1.68 

78.852 

61.853 

45.436 

95.218 

L21 

2.0 

1.6 

0.8 

1.68 

47.310 

40.701 

35.067 

159.857 

EL56 

5/6 

1.33 

1.6 

1.68 

56.741 

64.934 

6.653 

855.259 

EL65 

6/5 

1.6 

1.33 

1.68 

78.852 

65.678 

10.307 

356.416 

ML56 

5/6 

1.33 

1.6 

1.68 

56.741 

92.528 

46.025 

93.201 

ML65 

6/5 

1.6 

1.33 

1.68 

78.852 

105.559 

35.390 

157.635 

Xll 

1.0 

3.2 

3.2 

3.36 

756.981 

571.840 

193.961 

5.225 

X12 

0.5 

1.6 

3.2 

3.36 

189.721 

157.458 

47.869 

85.789 

X21 

2.0 

3.2 

1.6 

3.36 

378.481 

286.662 

65.393 

45.969 


Calculated at t = 1.4tci,h- 



20 



4 4.5 5 5.5 6 4 4.5 5 5.5 6 

log(T) log(T) 

Fig. 1.— Left panel: Cooling rate coefficient, L(T) = A/n'jj, for the standard cooling 
model (NEQmO, solid line), the collisional ionization equilibrium cooling model for a zero- 
metalicity gas (CIEmO, dotted line), and the non-equilibrium cooling model for a gas with 
the metalicity, Z = 2’q/IO (NEQml, dashed line) (see text). Right panel: Cooling time, 
tcooi = e/A, for a gas cooling from = 1.7 x 10®K with Uh = O.lcm”^ initially under the 
isobaric condition. The same line type is used for different cooling models as in the cooling 
rate coefficient. 
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[R^= 1Spherical Symmetric] 
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Fig. 2.— Time evolution of the gas density, p, the velocities, wr, and the pressure, P, at 
VPi,h = 0.2 (solid), 0.6 (dotted), 1.0 (dashed), 1.4 (long dashed) and 1.8 (dot and dashed line) 
for the spherical cloud with = 1.6/cooi (the Lll model). Distributions along the diagonal 
line of R = z in 2D box are shown, as a function of the radial distance r = vP^+^- 
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Fig. 3.— Time evolution of the density distribution in the prolate cloud with S = 5/6 and 
-Rc,max = 1-6/cool (the model L56). The contour levels for (p/po) = 2* (/ = 0,1, ••■10) are 
shown. Full images are generated from the simulation data covering one quadrant of the 
whole space by using the reflection symmetry along the z = 0 plane and along the ; 2 -axis. 
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[vfield : L56(Prolate) vs L65(Oblate)] 
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Fig. 5.— The velocity field overlapped on the density contour map for the prolate model 
L56 (left panels) and for the oblate model L65 (right panels) at t = 0.8fci,h (top panels) and 
at t = 1.6fci,h (bottom panels). 






















































































































































































































Fig. 8.— Density contour maps at t = 1.6tci,h for all the models with the standard cooling. 
The contour levels for (p/po) = 2* (i = 0,1, • ■ ■ 10) are shown. Each map is labeled with the 
model name (see Table 1). 
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Fig. 9.— Density line profiles of the cooled clonds for the S12 (top-left), S21 (top-right), 
M12 (npper middle-left) and M21 (npper middle-right) models at t = 1.2tci,h, and for the 
L12 (lower middle-left), L21 (lower middle-right), X12 (bottom-left), and X21 (bottom-right) 
models at t = 1.4tci,h- The dotted and dashed lines show the prohles along the R and 2 :-axes, 
respectively. The results of ID spherical simulations (the Sll, Mil, Lll, and Xll models) 
are shown with the solid lines for comparison. 
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Fig. 11.— Left panels: Mean density of the gas with T = 10"^K, < >, for prolate, spherical 

and oblate clouds (from bottom to top) shown as a function of time. The solid lines with 
circles are for the small cloud models, the dotted lines with triangles for the medium size 
cloud models, the dashed lines with squares for the large cloud models, and the long-dashed 
lines with crosses for the very large cloud models. Right panels: Total mass of the gas with 
T = 10"^K, M 4 , shown with the same line types and symbols as for < p 4 >. The spherical 
Jeans mass estimated with Tc = 10"^K and < >, Mj, is also shown in right panels with 

the same line types but without symbols. 
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